Module MCMC_data
  use netcdf

  Implicit none

  !********************* INPUTS [dimensions] ********************
  ! Dimensions:
  ! STATESIZE             = # of elements in emissions state vector, x_ap
  ! NMEASURE              = # of elements in measurement vector, z
  ! DIM1                  = maximum # of measurements that correspond to an element of the temporal model-measurement uncertainty (e.g., 100 if 100 observations out of nmeasuretotal correspond to element 1, and 80 correspond to element 2, etc)
  ! DIM2                  = # of elements of temporal model-measurement uncertainty estimated in the inversion (e.g., daytime and nighttime uncertainties over the month = 2 elements)
  ! DIM3                  = # of correlation parameters (e.g., tau or nu or rho) to be estimated for each of the spatial and temporal matrices (currently only value allowed is 1)
  ! NMEASURETOTAL         = maximum possible number of measurements during the time period of interest (e.g., 372 measurements in July if estimated every 2 hours).
  ! NUMSITES              = number of measurement sites in the inversion
  !                         NMEASUREMAX = NMEASURETOTAL*NUMSITES is computed in the code
  ! DIM4                  = number of gridcells in inversion domain

  ! Variables:
  ! NIT                   = Number of iterations [scalar]
  ! BURN_IN               = Length of burn-in chain [scalar]
  ! z                     = Vector of observations [nmeasure]
  ! H                     = sensitvity matrix [nmeasuremax x statesize]
  ! D                     = measurement uncertainty (nugget) matrix, assuming uncorrelated observations [nmeasure x nmeasure]

  ! X_AP                  = Prior parameter estimate vector [statesize]
  ! PDF_PARAM1            = Prior of primary parameter describing PDF (e.g. Mean) [statesize]
  ! PDF_PARAM2            = Prior of secondary parameter describing PDF (e.g. Standard Deviation) [statesize]
  ! SIGMA_Y               = Prior of model-measurement uncertainty in time [dim2]
  ! SIGMA_YS              = Prior of model-measurement uncertainty in space [currently fixed at numsites]
  ! TAU                   = Prior of autocorrelation timescale [dim3]
  ! NU                    = Prior of Matern shape parameter for spatial autocorrelation [dim3]
  ! RHO                   = Prior of Matern rate parameter for spatial autocorrelation [dim3]

  ! STEPSIZE              = Size of proposal distriubtion for x [statesize]
  ! STEPSIZE_PDF_PARAM1   = Size of proposal distriubtion for pdf_param1 [statesize]
  ! STEPSIZE_PDF_PARAM2   = Size of proposal distriubtion for pdf_param2 [statesize]
  ! STEPSIZE_SIGMA_Y      = Size of proposal distribution for sigma_y [dim2]
  ! STEPSIZE_SIGMA_YS     = Size of proposal distribution for sigma_ys [numsites]
  ! STEPSIZE_TAU          = Size of proposal distribution for tau [dim3]
  ! STEPSIZE_NU           = Size of proposal distribution for nu [dim3]
  ! STEPSIZE_RHO          = Size of proposal distribution for rho [dim3]
  ! STEPSIZE_Y            = Size of proposal distribution for y [nmeasuremax]

  ! PDF type parameters ( 0 = LOGNORMAL, 1 = EXPONENTIAL, 2 = GAUSSIAN, 3 = UNIFORM) [statesize]

  ! X_PDF                 = Type of PDF associated with each element of x[statesize]
  ! PDF_PARAM1_PDF        = Type of PDF associated with each element of pdf_param1 [statesize]
  ! PDF_PARAM2_PDF        = Type of PDF associated with each element of pdf_param2 [statesize]
  ! SIGMA_Y_PDF           = Type of PDF associated with each element of sigma_y [dim2]
  ! SIGMA_YS_PDF          = Type of PDF associated with each element of sigma_ys [numsites]
  ! TAU_PDF               = Type of PDF associated with each element of tau [dim3]
  ! NU_PDF                = Type of PDF associated with each element of nu [dim3]
  ! RHO_PDF               = Type of PDF associated with each element of rho [dim3]

  ! PDF_HYPERPARAM1       = Parameter governing PDF of pdf_param1 (e.g., standard deviation in pdf_param1 "hyper-PDF") [statesize]
  ! PDF_HYPERPARAM2       = Parameter governing PDF of pdf_param2 (e.g., uncertainty in pdf_param2 "hyper-PDF") [statesize]
  ! SIGMA_Y_HYPERPARAM    = Parameter governing PDF of sigma_y (e.g., uncertainty in sigma_y "hyper-PDF") [dim2]
  ! SIGMA_YS_HYPERPARAM   = Parameter governing PDF of sigma_ys (e.g., uncertainty in sigma_ys "hyper-PDF") [numsites]
  ! TAU_HYPERPARAM1       = Parameter governing PDF of tau [dim3]
  ! TAU_HYPERPARAM2       = Parameter governing PDF of tau [dim3]
  ! NU_HYPERPARAM1        = Parameter1 governing PDF of nu [dim3]
  ! RHO_HYPERPARAM1       = Parameter1 governing PDF of rho [dim3]
  ! NU_HYPERPARAM2        = Parameter2 governing PDF of nu [dim3]
  ! RHO_HYPERPARAM2       = Parameter2 governing PDF of rho [dim3]

  ! T_INDICES             = Array containing indices of measurements that correspond to each element of sigma_y (e.g., 1st column corresponds to the indices governed by first element of sigma_y, 2nd column  corresponds to indices governed by the second element of sigma_y, and so forth) [dim1 x dim2]

  ! DELTATIME             = Array of 'delta_t' in same units as TAU to calculate temporal autocorrelation (e.g., if tau is in days, measurements that occur 3 hours apart would have a value of 0.125 days, 6 hours apart would have a value of 0.5, etc.) [nmeasuretotal x nmeasuretotal]

  ! DISTANCE              = Array of Euclidean distance between sites in same units as rho to calculate spatial autocorrelation  [numsites x numsites]

  ! DATENUMBER            = Vector of serial date number (or equivalent) for the time period (e.g., July 2012 estimated every 3 hours would have 248 datenumbers). This is not used in the MCMC routine but is saved into the outputs for the convenience of analysis. Ignore if desired.  [nmeasuretotal]

  ! TIMEINDEX_NONZERO     = Vector of indices corresponding to when a measurement was avaiable during time period (1:nmeasuremax).  [nmeasure]

  ! N_OBS                 = Vector of number observations per site (not used in analysis)  [numsites]

  ! KRON_FLAG             = Flag for orientation of y,H and kronecker product. Value = 0 if spatially stacked, = 1 if temporally stacked. [Scalar]

  ! NUMTHREADS            = Number of threads to use for multithreaded applications

  ! SITESPRESENT          = Index corresponding to which sites are present using whatever numbering system the user has implemented (not used in analysis but important for tracking which sites are used). For example if Mace Head = 1, Ridge Hill = 2, Tacolneston = 3, Angus = 4 then if only MHD and TAC are available for the analysis then this would be [1 3] [numsites]

  ! DISTRIBUTION          = Percentage of emissions from grid cell contained in larger region that has been aggregated for the inversion (not used in the analysis) [dim4]

  !! *************** OUTPUTS [dimensions] ***************************
  ! X_IT                  = Array containing value of x at each iteration (nIt x statesize)
  ! PDF_PARAM1_IT         = Array containing value of pdf_param1 at each iteration (nIt x statesize)
  ! PDF_PARAM2_IT         = Array containing value of pdf_param2 at each iteration (nIt x statesize)
  ! SIGMA_Y_IT            = Array containing value of sigma_y at each iteration (nIt x dim2)
  ! SIGMA_YS_IT           = Array containing value of sigma_ys at each iteration (nIt x numsites)
  ! TAU_IT                = Array containing value of tau at each iteration (nIt x dim3)
  ! NU_IT                 = Array containing value of nu at each iteration (nIt x dim3)
  ! RHO_IT                = Array containing value of rho at each iteration (nIt x dim3)
  ! Y_IT                  = Array containing value of y at each iteration (nIt x nmeasuremax)
  ! ACCEPTANCE            = Vector of acceptance ratio for each element of x/pdf_param1/pdf_param2 [statesize]
  ! ACCEPTANCE_SIGMAY     = Vector of acceptance ratio for each element of sigma_y [dim2]
  ! ACCEPTANCE_SIGMAYS    = Vector of acceptance ratio for each element of sigma_ys [numsites]
  ! ACCEPTANCE_TAU        = Vector of acceptance ratio for each element of tau [dim3]
  ! ACCEPTANCE_NU         = Vector of acceptance ratio for each element of nu and rho [dim3]
  ! ACCEPTANCE_Y          = Vector of acceptance ratio for each element of y [nmeasuremax]
  ! H                     = Same as input - restored for convenience of analyzing output
  ! Z                     = Same as input - restored for convenience of analyzing output
  ! DATENUMBER            = Same as input - restored for convenience of analyzing output
  ! TIMEINDEX_NONZERO     = Same as input - restored for convenience of analyzing output
  ! X_AP                  = Same as input - restored for convenience of analyzing output
  ! N_OBS                 = Same as input - restored for convenience of analyzing output
  ! sitespresent          = Same as input - restored for convenience of analyzing output
  ! deltatime             = Same as input - restored for convenience of analyzing output
  ! distance              = Same as input - restored for convenience of analyzing output
  ! T_indices             = Same as input - restored for convenience of analyzing output
  ! DISTRIBUTION          = Same as input - restored for convenience of analyzing output
  ! D                     = Same as input - restored for convenience of analyzing output

  ! Variable declarations

  ! dimensions
  integer                           ::  nIt, statesize, nmeasure, nmeasuretotal, nmeasuremax, numsites, burn_in, dim1, dim2, &
       dim3, kron_flag, numthreads, dim4

  ! counters
  integer                           ::  status, tstart, tend, &
       reject, reject_sigma_y, reject_tau , reject_nu, reject_y, reject_sigma_ys, omp_get_num_threads
  ! input and output variable and dimension IDs
  integer                           ::  inputID, outputID, burn_inID, nItID, statesizeID, nMeasureID, x_apID, x_pdfID, &
       pdf_param1ID, pdf_param2ID, stepsizeID, zID, HID, x_itID, y_itID, acceptID, &
       statesizeDID, nMeasureDID, nItDID, stepsizepdfparam1ID, stepsizepdfparam2ID, &
       pdf_hyperparam1ID, pdf_hyperparam2ID, pdf_param1_pdfID, pdf_param2_pdfID, &
       sigmayID, TindicesID, sigmayhyperparamID,  stepsizesigmayID, sigmaypdfID, &
       dim1ID, dim2ID, dim3ID, pdf_param1_itID, pdf_param2_itID, sigma_y_itID, &
       dim1DID, dim2DID, dim3DID, acceptsigmayID, acceptyID, deltatimeID, distanceID, &
       tauID, taupdfID, tau_itID, accepttauID, stepsizetauID, tauhyperparam1ID, &
       tauhyperparam2ID, HoID, zoID, nmeasuretotalID, datenumberID, timeindexnonzeroID, &
       datenumberoID, timeindexnonzerooID, stepsizenuID, stepsizerhoID, nuID, rhoID, &
       nuhyperparam1ID, rhohyperparam1ID, nuhyperparam2ID, rhohyperparam2ID, nupdfID, &
       rhopdfID, acceptnuID, nu_itID, rho_itID, numsitesID, DID, stepsizeyID, &
       nmeasuremaxDID, xapoID, nobsID, nobsoID, numsitesDID, sigmaysID, &
       sigmays_itID, acceptsigmaysID, stepsizesigmaysID, sigmayshyperparamID, &
       sigmayspdfID, kron_flagID, numthreadsID, sitespresentID, sitespresentoID, &
       deltatimeoID, distanceoID, TindicesoID, dim4ID, distributionID, dim4DID, &
       distributionoID, DoID, nmeasuretotalDID

  integer                           ::  tau_pdf, nu_pdf, rho_pdf, count_rate, count_0, count_1, count_max

  integer ,allocatable, dimension(:) ::  x_pdf, pdf_param1_pdf, pdf_param2_pdf, sigma_y_pdf, sigma_ys_pdf, &
       reject_vector, reject_sigma_y_vector, reject_sigma_ys_vector, reject_y_vector, &
       timeindex_nonzero, nobs, sitespresent, num_T_indices, timeindex

  integer, allocatable, dimension(:,:) :: T_indices

  real                              ::  n0T, n1T, m0T, m1T, &
       p0_tau, p0_nu, p0_rho, p1, &
       tau, tau_current, stepsize_tau, tau_hyperparam1, &
       tau_hyperparam2, accept_tau, accept_nu, nu, nu_current, stepsize_nu, &
       nu_hyperparam1, nu_hyperparam2, rho, rho_current, stepsize_rho, &
       rho_hyperparam1, rho_hyperparam2, detval_current,detval_new, detval, detval2, &
       detval_T_current, detval_T_new, detval_S_current, detval_S_new, ga, aa, bb, &
       t0

  real,allocatable, dimension(:)    ::  stepsize, stepsize_pdf_param1, stepsize_pdf_param2, stepsize_sigma_y, &
       stepsize_sigma_ys, stepsize_y, p0_x,p0_pdf_param1, p0_pdf_param2, &
       p0_sigma_y, p0_sigma_ys, x, x_ap, z, y, dy, y_small, pdf_param1, pdf_param2, &
       pdf_param1_current, pdf_param1_new, pdf_param2_current, pdf_param2_new, &
       pdf_hyperparam1, pdf_hyperparam2, y_current, y_new, y_error_t, dumx, dumy, dum2, &
       accept_vector, accept_sigma_y_vector, accept_sigma_ys_vector, accept_y_vector, &
       sigma_y, sigma_y_current, sigma_y_hyperparam, autocorr_vec, &
       sigma_ys, sigma_ys_current, sigma_ys_hyperparam, &
       datenumber, distribution, n0, n1, m0, m1, C, C2

  real,allocatable, dimension(:,:)  ::  H, D, Dinv, deltatime, distance, x_it, pdf_param1_it, pdf_param2_it, sigma_y_it, &
       tau_it, nu_it, rho_it, y_it, sigma_ys_it,  Rinv_current, &
       T_current, T_new, S_current, S_new, Sinv_current, Sinv_new, Tinv_current, &
       Tinv_new, dum3

  real (kind =8) nu_double, rho_double, alpha, ga_double, arg_dum, k_arg_dum

  real (kind =8), allocatable, dimension(:,:)   :: distance_double, arg

  integer ( kind = 4 ) nb, ize, ncalc


Contains

  subroutine allocate_mem

    integer :: errstat
    integer, parameter :: r=4, r8=8, i=4
    integer :: bytes

    bytes = 0

    allocate(z(nMeasure),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate z'
       stop
    endif
    bytes = bytes + nMeasure * r

    allocate(y(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate y'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(nobs(numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate nobs'
       stop
    endif
    bytes = bytes + numsites * i

    allocate(sitespresent(numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate sitespresent'
       stop
    endif
    bytes = bytes + numsites * i

    allocate(y_current(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate y_current'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(y_new(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate y_new'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(y_small(nmeasure),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate y_small'
       stop
    endif
    bytes = bytes + nmeasure * r

    allocate(n0(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate n0'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(n1(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate n1'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(m0(nmeasure),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate m0'
       stop
    endif
    bytes = bytes + nmeasure * r

    allocate(m1(nmeasure),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate m1'
       stop
    endif
    bytes = bytes + nmeasure * r

    allocate(dy(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate dy'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(H(nmeasuremax,statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate H'
       stop
    endif
    bytes = bytes + nmeasuremax * statesize * r

    allocate(D(nMeasure,nMeasure),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate D'
       stop
    endif
    bytes = bytes + nMeasure * nMeasure * r

    allocate(Dinv(nMeasure,nMeasure),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate Dinv'
       stop
    endif
    bytes = bytes + nMeasure * nMeasure * r

    allocate(T_indices(dim1,dim2),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate T_indices'
       stop
    endif
    bytes = bytes + dim1 * dim2 * i

    allocate(deltatime(nmeasuretotal,nmeasuretotal),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate deltatime'
       stop
    endif
    bytes = bytes + nmeasuretotal * nmeasuretotal * r

    allocate(distance(numsites,numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate distance'
       stop
    endif
    bytes = bytes + numsites * numsites * r

    allocate(distance_double(numsites,numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate distance_double'
       stop
    endif
    bytes = bytes + numsites * numsites * r8

    allocate(arg(numsites,numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate arg'
       stop
    endif
    bytes = bytes + numsites * numsites * r8

    allocate(x_ap(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate x_ap'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(x(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate x'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(pdf_param1(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_param1'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(pdf_param1_current(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_param1_current'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(pdf_param1_new(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_param1_new'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(pdf_param2(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_param2'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(pdf_param2_current(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_param2_current'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(pdf_param2_new(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_param1_new'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(sigma_y(dim2),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate sigma_y'
       stop
    endif
    bytes = bytes + dim2 * r

    allocate(sigma_y_current(dim2),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate sigma_y_current'
       stop
    endif
    bytes = bytes + dim2 * r

    allocate(sigma_ys(numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate sigma_ys'
       stop
    endif
    bytes = bytes + numsites * r

    allocate(sigma_ys_current(numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate sigma_ys_current'
       stop
    endif
    bytes = bytes + numsites * r

    allocate(stepsize(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate stepsize'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(stepsize_pdf_param1(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate stepsize_pdf_param1'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(stepsize_pdf_param2(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate stepsize_pdf_param2'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(stepsize_sigma_y(dim2),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate stepsize_sigma_y'
       stop
    endif
    bytes = bytes + dim2 * r

    allocate(stepsize_sigma_ys(numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate stepsize_sigma_ys'
       stop
    endif
    bytes = bytes + numsites * r

    allocate(stepsize_y(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate stepsize_y'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(x_pdf(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate x_pdf'
       stop
    endif
    bytes = bytes + statesize * i

    allocate(pdf_param1_pdf(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_param1_pdf'
       stop
    endif
    bytes = bytes + statesize * i

    allocate(pdf_param2_pdf(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_param2_pdf'
       stop
    endif
    bytes = bytes + statesize * i

    allocate(sigma_y_pdf(dim2),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate sigma_y_pdf'
       stop
    endif
    bytes = bytes + dim2 * i

    allocate(sigma_ys_pdf(numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate sigma_ys_pdf'
       stop
    endif
    bytes = bytes + numsites * i

    allocate(pdf_hyperparam1(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_hyperparam1'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(pdf_hyperparam2(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_hyperparam2'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(sigma_y_hyperparam(dim2),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate sigma_y_hyperparam1'
       stop
    endif
    bytes = bytes + dim2 * r

    allocate(sigma_ys_hyperparam(numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate sigma_ys_hyperparam1'
       stop
    endif
    bytes = bytes + numsites * r

    allocate(p0_x(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate p0_x'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(p0_pdf_param1(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate p0_pdf_param1'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(p0_pdf_param2(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate p0_pdf_param2'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(p0_sigma_y(dim2),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate p0_sigma_y'
       stop
    endif
    bytes = bytes + dim2 * r

    allocate(p0_sigma_ys(numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate p0_sigma_ys'
       stop
    endif
    bytes = bytes + numsites * r

    allocate(y_error_t(nmeasuretotal),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate y_error_t array'
       stop
    endif
    bytes = bytes + nmeasuretotal * r

    allocate(dumx(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate dumx'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(dumy(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate dumy'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(dum2(nmeasure),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate dum2'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(C(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate C'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(C2(nmeasure),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate C2'
       stop
    endif
    bytes = bytes + nmeasure * r

    allocate(T_current(nmeasuretotal,nmeasuretotal),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate T_current'
       stop
    endif
    bytes = bytes + nmeasuretotal * nmeasuretotal * r

    allocate(T_new(nmeasuretotal,nmeasuretotal),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate T_new'
       stop
    endif
    bytes = bytes + nmeasuretotal * nmeasuretotal * r

    allocate(Tinv_current(nmeasuretotal,nmeasuretotal),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate Tinv_current'
       stop
    endif
    bytes = bytes + nmeasuretotal * nmeasuretotal * r

    allocate(Tinv_new(nmeasuretotal,nmeasuretotal),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate Tinv_new'
       stop
    endif
    bytes = bytes + nmeasuretotal * nmeasuretotal * r

    allocate(S_current(numsites,numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate S_current'
       stop
    endif
    bytes = bytes + numsites * numsites * r

    allocate(S_new(numsites,numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate S_new'
       stop
    endif
    bytes = bytes + numsites * numsites * r

    allocate(Sinv_current(numsites,numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate Sinv_current'
       stop
    endif
    bytes = bytes + numsites * numsites * r

    allocate(Sinv_new(numsites,numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate Sinv_new'
       stop
    endif
    bytes = bytes + numsites * numsites * r

    allocate(Rinv_current(nmeasuremax,nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate Rinv_current'
       stop
    endif
    bytes = bytes + nmeasuremax * nmeasuremax * r

    allocate(autocorr_vec(nmeasuretotal),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate autocorr_vec'
       stop
    endif
    bytes = bytes + nmeasuretotal * r

    allocate(x_it(nIt,statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate x_it'
       stop
    endif
    bytes = bytes + nIt * statesize * r

    allocate(pdf_param1_it(nIt,statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_param1_it'
       stop
    endif
    bytes = bytes + nIt * statesize * r

    allocate(pdf_param2_it(nIt,statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate pdf_param2_it'
       stop
    endif
    bytes = bytes + nIt * statesize * r

    allocate(sigma_y_it(nIt,dim2),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate sigma_y_it'
       stop
    endif
    bytes = bytes + nIt * dim2 * r

    allocate(tau_it(nIT,dim3),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate tau_it'
       stop
    endif
    bytes = bytes + nIT * dim3 * r

    allocate(nu_it(nIT,dim3),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate nu_it'
       stop
    endif
    bytes = bytes + nIT * dim3 * r

    allocate(rho_it(nIT,dim3),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate rho_it'
       stop
    endif
    bytes = bytes + nIT * dim3 * r

    allocate(sigma_ys_it(nIT,numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate sigma_ys_it'
       stop
    endif
    bytes = bytes + nIT * numsites * r

    allocate(y_it(nIT,nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate y_it'
       stop
    endif
    bytes = bytes + nIT * nmeasuremax * r

    allocate(reject_vector(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate reject_vector'
       stop
    endif
    bytes = bytes + statesize * i

    allocate(reject_sigma_y_vector(dim2),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate reject_sigma_y_vector'
       stop
    endif
    bytes = bytes + dim2 * i

    allocate(reject_sigma_ys_vector(numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate reject_sigma_ys_vector'
       stop
    endif
    bytes = bytes + numsites * i

    allocate(reject_y_vector(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate reject_y_vector'
       stop
    endif
    bytes = bytes + nmeasuremax * i

    allocate(accept_vector(statesize),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate accept_vector'
       stop
    endif
    bytes = bytes + statesize * r

    allocate(accept_sigma_y_vector(dim2),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate accept_sigma_y_vector'
       stop
    endif
    bytes = bytes + dim2 * r

    allocate(accept_sigma_ys_vector(numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate accept_sigma_ys_vector'
       stop
    endif
    bytes = bytes + numsites * r

    allocate(accept_y_vector(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate accept_y_vector'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(datenumber(nmeasuretotal),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate datenumber'
       stop
    endif
    bytes = bytes + nmeasuretotal * r

    allocate(timeindex(nmeasuremax),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate timeindex'
       stop
    endif
    bytes = bytes + nmeasuremax * r

    allocate(timeindex_nonzero(nmeasure),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate timeindex_nonzero'
       stop
    endif
    bytes = bytes + nmeasure * r

    allocate(distribution(dim4),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate distribution'
       stop
    endif
    bytes = bytes + dim4 * r

    allocate(dum3(nmeasuretotal,numsites),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate dum3'
       stop
    endif
    bytes = bytes + nmeasuretotal * numsites * r

    allocate(num_T_indices(dim2),stat=errstat)
    if (errstat /= 0) then
       write (*,*) 'ERROR: could not allocate num_T_indices'
       stop
    endif
    bytes = bytes + dim2 * i

    write(*,*) "Allocated :",bytes/1024/1024,"MB"

  end subroutine allocate_mem

  subroutine open_files

    STATUS = NF90_OPEN('./inputs/input_file.nc', NF90_NOWRITE, ncid=inputID)
    STATUS = NF90_CREATE ('./outputs/output_file.nc', cmode=0, ncid=outputID)

  end subroutine open_files

  subroutine read_params

    STATUS = NF90_INQ_VARID (inputID, 'burn_in', burn_inID)
    STATUS = NF90_INQ_VARID (inputID, 'nit', nItID)
    STATUS = NF90_INQ_VARID (inputID, 'statesize', statesizeID)
    STATUS = NF90_INQ_VARID (inputID, 'nmeasure', nMeasureID)
    STATUS = NF90_INQ_VARID (inputID, 'dim1', dim1ID)
    STATUS = NF90_INQ_VARID (inputID, 'dim2', dim2ID)
    STATUS = NF90_INQ_VARID (inputID, 'dim3', dim3ID)
    STATUS = NF90_INQ_VARID (inputID, 'nmeasuretotal', nmeasuretotalID)
    STATUS = NF90_INQ_VARID (inputID, 'numsites', numsitesID)
    STATUS = NF90_INQ_VARID (inputID, 'dim4', dim4ID)

    STATUS = NF90_GET_VAR (inputID, burn_inID, burn_in)
    STATUS = NF90_GET_VAR (inputID, nItID, nIt)
    STATUS = NF90_GET_VAR (inputID, statesizeID, statesize)
    STATUS = NF90_GET_VAR (inputID, nMeasureID, nMeasure)
    STATUS = NF90_GET_VAR (inputID, dim1ID, dim1)
    STATUS = NF90_GET_VAR (inputID, dim2ID, dim2)
    STATUS = NF90_GET_VAR (inputID, dim3ID, dim3)
    STATUS = NF90_GET_VAR (inputID, nmeasuretotalID, nmeasuretotal)
    STATUS = NF90_GET_VAR (inputID, numsitesID, numsites)
    STATUS = NF90_GET_VAR (inputID, dim4ID, dim4)

    nmeasuremax = nmeasuretotal*numsites

  end subroutine read_params

  subroutine read_data

    STATUS = NF90_INQ_VARID (inputID, 'x_ap', x_apID)
    STATUS = NF90_INQ_VARID (inputID, 'stepsize', stepsizeID)
    STATUS = NF90_INQ_VARID (inputID, 'stepsize_pdf_param1', stepsizepdfparam1ID)
    STATUS = NF90_INQ_VARID (inputID, 'stepsize_pdf_param2', stepsizepdfparam2ID)
    STATUS = NF90_INQ_VARID (inputID, 'stepsize_sigma_y', stepsizesigmayID)
    STATUS = NF90_INQ_VARID (inputID, 'stepsize_tau', stepsizetauID)
    STATUS = NF90_INQ_VARID (inputID, 'stepsize_nu', stepsizenuID)
    STATUS = NF90_INQ_VARID (inputID, 'stepsize_rho', stepsizerhoID)
    STATUS = NF90_INQ_VARID (inputID, 'stepsize_y', stepsizeyID)
    STATUS = NF90_INQ_VARID (inputID, 'stepsize_sigma_ys', stepsizesigmaysID)
    STATUS = NF90_INQ_VARID (inputID, 'pdf_param1', pdf_param1ID)
    STATUS = NF90_INQ_VARID (inputID, 'pdf_param2', pdf_param2ID)
    STATUS = NF90_INQ_VARID (inputID, 'sigma_y', sigmayID)
    STATUS = NF90_INQ_VARID (inputID, 'sigma_ys', sigmaysID)
    STATUS = NF90_INQ_VARID (inputID, 'tau', tauID)
    STATUS = NF90_INQ_VARID (inputID, 'nu', nuID)
    STATUS = NF90_INQ_VARID (inputID, 'rho', rhoID)
    STATUS = NF90_INQ_VARID (inputID, 'pdf_hyperparam1', pdf_hyperparam1ID)
    STATUS = NF90_INQ_VARID (inputID, 'pdf_hyperparam2', pdf_hyperparam2ID)
    STATUS = NF90_INQ_VARID (inputID, 'sigma_y_hyperparam', sigmayhyperparamID)
    STATUS = NF90_INQ_VARID (inputID, 'tau_hyperparam1', tauhyperparam1ID)
    STATUS = NF90_INQ_VARID (inputID, 'tau_hyperparam2', tauhyperparam2ID)
    STATUS = NF90_INQ_VARID (inputID, 'sigma_ys_hyperparam', sigmayshyperparamID)
    STATUS = NF90_INQ_VARID (inputID, 'nu_hyperparam1', nuhyperparam1ID)
    STATUS = NF90_INQ_VARID (inputID, 'rho_hyperparam1', rhohyperparam1ID)
    STATUS = NF90_INQ_VARID (inputID, 'nu_hyperparam2', nuhyperparam2ID)
    STATUS = NF90_INQ_VARID (inputID, 'rho_hyperparam2', rhohyperparam2ID)
    STATUS = NF90_INQ_VARID (inputID, 'z', zID)
    STATUS = NF90_INQ_VARID (inputID, 'nobs', nobsID)
    STATUS = NF90_INQ_VARID (inputID, 'H', HID)
    STATUS = NF90_INQ_VARID (inputID, 'D', DID)
    STATUS = NF90_INQ_VARID (inputID, 'x_pdf', x_pdfID)
    STATUS = NF90_INQ_VARID (inputID, 'pdf_param1_pdf', pdf_param1_pdfID)
    STATUS = NF90_INQ_VARID (inputID, 'pdf_param2_pdf', pdf_param2_pdfID)
    STATUS = NF90_INQ_VARID (inputID, 'tau_pdf', taupdfID)
    STATUS = NF90_INQ_VARID (inputID, 'nu_pdf', nupdfID)
    STATUS = NF90_INQ_VARID (inputID, 'rho_pdf', rhopdfID)
    STATUS = NF90_INQ_VARID (inputID, 'sigma_y_pdf', sigmaypdfID)
    STATUS = NF90_INQ_VARID (inputID, 'sigma_ys_pdf', sigmayspdfID)
    STATUS = NF90_INQ_VARID (inputID, 'T_indices', TindicesID)
    STATUS = NF90_INQ_VARID (inputID, 'deltatime', deltatimeID)
    STATUS = NF90_INQ_VARID (inputID, 'distance', distanceID)
    STATUS = NF90_INQ_VARID (inputID, 'datenumber', datenumberID)
    STATUS = NF90_INQ_VARID (inputID, 'timeindex_nonzero', timeindexnonzeroID)
    STATUS = NF90_INQ_VARID (inputID, 'kron_flag', kron_flagID)
    STATUS = NF90_INQ_VARID (inputID, 'numthreads', numthreadsID)
    STATUS = NF90_INQ_VARID (inputID, 'sitespresent', sitespresentID)
    STATUS = NF90_INQ_VARID (inputID, 'distribution', distributionID)

    STATUS = NF90_GET_VAR (inputID, x_apID, x_ap)
    STATUS = NF90_GET_VAR (inputID, stepsizeID, stepsize)
    STATUS = NF90_GET_VAR (inputID, stepsizepdfparam1ID, stepsize_pdf_param1)
    STATUS = NF90_GET_VAR (inputID, stepsizepdfparam2ID, stepsize_pdf_param2)
    STATUS = NF90_GET_VAR (inputID, stepsizesigmayID, stepsize_sigma_y)
    STATUS = NF90_GET_VAR (inputID, stepsizetauID, stepsize_tau)
    STATUS = NF90_GET_VAR (inputID, stepsizenuID, stepsize_nu)
    STATUS = NF90_GET_VAR (inputID, stepsizerhoID, stepsize_rho)
    STATUS = NF90_GET_VAR (inputID, stepsizeyID, stepsize_y)
    STATUS = NF90_GET_VAR (inputID, stepsizesigmaysID, stepsize_sigma_ys)
    STATUS = NF90_GET_VAR (inputID, pdf_param1ID, pdf_param1)
    STATUS = NF90_GET_VAR (inputID, pdf_param2ID, pdf_param2)
    STATUS = NF90_GET_VAR (inputID, sigmayID, sigma_y)
    STATUS = NF90_GET_VAR (inputID, sigmaysID, sigma_ys)
    STATUS = NF90_GET_VAR (inputID, tauID, tau)
    STATUS = NF90_GET_VAR (inputID, nuID, nu)
    STATUS = NF90_GET_VAR (inputID, rhoID, rho)
    STATUS = NF90_GET_VAR (inputID, pdf_hyperparam1ID, pdf_hyperparam1)
    STATUS = NF90_GET_VAR (inputID, pdf_hyperparam2ID, pdf_hyperparam2)
    STATUS = NF90_GET_VAR (inputID, sigmayhyperparamID, sigma_y_hyperparam)
    STATUS = NF90_GET_VAR (inputID, tauhyperparam1ID, tau_hyperparam1)
    STATUS = NF90_GET_VAR (inputID, tauhyperparam2ID, tau_hyperparam2)
    STATUS = NF90_GET_VAR (inputID, nuhyperparam1ID, nu_hyperparam1)
    STATUS = NF90_GET_VAR (inputID, rhohyperparam1ID, rho_hyperparam1)
    STATUS = NF90_GET_VAR (inputID, nuhyperparam2ID, nu_hyperparam2)
    STATUS = NF90_GET_VAR (inputID, rhohyperparam2ID, rho_hyperparam2)
    STATUS = NF90_GET_VAR (inputID, sigmayshyperparamID, sigma_ys_hyperparam)
    STATUS = NF90_GET_VAR (inputID, zID, z)
    STATUS = NF90_GET_VAR (inputID, nobsID, nobs)
    STATUS = NF90_GET_VAR (inputID, HID, H)
    STATUS = NF90_GET_VAR (inputID, DID, D)
    STATUS = NF90_GET_VAR (inputID, x_pdfID, x_pdf)
    STATUS = NF90_GET_VAR (inputID, pdf_param1_pdfID, pdf_param1_pdf)
    STATUS = NF90_GET_VAR (inputID, pdf_param2_pdfID, pdf_param2_pdf)
    STATUS = NF90_GET_VAR (inputID, taupdfID, tau_pdf)
    STATUS = NF90_GET_VAR (inputID, nupdfID, nu_pdf)
    STATUS = NF90_GET_VAR (inputID, rhopdfID, rho_pdf)
    STATUS = NF90_GET_VAR (inputID, sigmaypdfID, sigma_y_pdf)
    STATUS = NF90_GET_VAR (inputID, sigmayspdfID, sigma_ys_pdf)
    STATUS = NF90_GET_VAR (inputID, TindicesID, T_indices)
    STATUS = NF90_GET_VAR (inputID, deltatimeID, deltatime)
    STATUS = NF90_GET_VAR (inputID, distanceID, distance)
    STATUS = NF90_GET_VAR (inputID, datenumberID, datenumber)
    STATUS = NF90_GET_VAR (inputID, timeindexnonzeroID, timeindex_nonzero)
    STATUS = NF90_GET_VAR (inputID, kron_flagID, kron_flag)
    STATUS = NF90_GET_VAR (inputID, numthreadsID, numthreads)
    STATUS = NF90_GET_VAR (inputID, sitespresentID, sitespresent)
    STATUS = NF90_GET_VAR (inputID, distributionID, distribution)

    STATUS = nf90_close(inputID)

  end subroutine read_data

  subroutine write_data

    STATUS = NF90_DEF_DIM (outputID, 'nit', nIt, nItDID)
    STATUS = NF90_DEF_DIM (outputID, 'statesize', statesize, statesizeDID)
    STATUS = NF90_DEF_DIM (outputID, 'nmeasure', nMeasure, nMeasureDID)
    STATUS = NF90_DEF_DIM (outputID, 'dim1', dim1, dim1DID)
    STATUS = NF90_DEF_DIM (outputID, 'dim2', dim2, dim2DID)
    STATUS = NF90_DEF_DIM (outputID, 'dim3', dim3, dim3DID)
    STATUS = NF90_DEF_DIM (outputID, 'nmeasuretotal', nmeasuretotal, nmeasuretotalDID)
    STATUS = NF90_DEF_DIM (outputID, 'nmeasuremax', nmeasuremax, nmeasuremaxDID)
    STATUS = NF90_DEF_DIM (outputID, 'numsites', numsites, numsitesDID)
    STATUS = NF90_DEF_DIM (outputID, 'dim4', dim4, dim4DID)


    STATUS = NF90_DEF_VAR (outputID, 'x_it', NF90_DOUBLE, (/ nItDID, statesizeDID/), x_itID)
    STATUS = NF90_DEF_VAR (outputID, 'pdf_param1_it', NF90_DOUBLE, (/ nItDID, statesizeDID/), pdf_param1_itID)
    STATUS = NF90_DEF_VAR (outputID, 'pdf_param2_it', NF90_DOUBLE, (/ nItDID, statesizeDID/), pdf_param2_itID)
    STATUS = NF90_DEF_VAR (outputID, 'sigma_y_it', NF90_DOUBLE, (/ nItDID, dim2DID/), sigma_y_itID)
    STATUS = NF90_DEF_VAR (outputID, 'sigma_ys_it', NF90_DOUBLE, (/ nItDID, numsitesDID/), sigmays_itID)
    STATUS = NF90_DEF_VAR (outputID, 'tau_it', NF90_DOUBLE, nItDID, tau_itID)
    STATUS = NF90_DEF_VAR (outputID, 'nu_it', NF90_DOUBLE, nItDID, nu_itID)
    STATUS = NF90_DEF_VAR (outputID, 'rho_it', NF90_DOUBLE, nItDID, rho_itID)
    STATUS = NF90_DEF_VAR (outputID, 'y_it', NF90_DOUBLE, (/ nItDID, nmeasuremaxDID/), y_itID)
    STATUS = NF90_DEF_VAR (outputID, 'acceptance', NF90_DOUBLE, statesizeDID, acceptID)
    STATUS = NF90_DEF_VAR (outputID, 'acceptance_sigmay', NF90_DOUBLE, dim2DID, acceptsigmayID)
    STATUS = NF90_DEF_VAR (outputID, 'acceptance_sigmays', NF90_DOUBLE, dim3DID , acceptsigmaysID)
    STATUS = NF90_DEF_VAR (outputID, 'acceptance_tau', NF90_DOUBLE, dim3DID , accepttauID)
    STATUS = NF90_DEF_VAR (outputID, 'acceptance_nu', NF90_DOUBLE, dim3DID , acceptnuID)
    STATUS = NF90_DEF_VAR (outputID, 'acceptance_y', NF90_DOUBLE, nmeasuremaxDID, acceptyID)
    STATUS = NF90_DEF_VAR (outputID, 'H', NF90_DOUBLE, (/ nMeasuremaxDID, statesizeDID/), HoID) ! include for convenience
    STATUS = NF90_DEF_VAR (outputID, 'z', NF90_DOUBLE, nMeasureDID, zoID) ! include for convenience
    STATUS = NF90_DEF_VAR (outputID, 'datenumber', NF90_DOUBLE, nmeasuretotalDID, datenumberoID) ! include for convenience
    STATUS = NF90_DEF_VAR (outputID, 'timeindex_nonzero', NF90_DOUBLE, nMeasureDID, timeindexnonzerooID) ! include for convenience
    STATUS = NF90_DEF_VAR (outputID, 'x_ap', NF90_DOUBLE, statesizeDID, xapoID) ! include for convenience
    STATUS = NF90_DEF_VAR (outputID, 'nobs', NF90_DOUBLE, numsitesDID, nobsoID) ! include for convenience
    STATUS = NF90_DEF_VAR (outputID, 'sitespresent', NF90_DOUBLE, numsitesDID, sitespresentoID) ! include for convenience
    STATUS = NF90_DEF_VAR (outputID, 'deltatime', NF90_DOUBLE, (/ nmeasuretotalDID, nmeasuretotalDID/), deltatimeoID) ! include for convenience
    STATUS = NF90_DEF_VAR (outputID, 'distance', NF90_DOUBLE, (/ numsitesDID, numsitesDID/), distanceoID) ! include for convenience
    STATUS = NF90_DEF_VAR (outputID, 'T_indices', NF90_DOUBLE, (/ dim1DID, dim2DID/), TindicesoID) ! include for convenience
    STATUS = NF90_DEF_VAR (outputID, 'distribution', NF90_DOUBLE, dim4DID, distributionoID) ! include for convenience
    STATUS = NF90_DEF_VAR (outputID, 'D', NF90_DOUBLE, (/ nMeasureDID, nMeasureDID/), DoID) ! include for convenience
    STATUS = NF90_ENDDEF (outputID)

    STATUS = NF90_PUT_VAR (outputID, x_itID, x_it)
    STATUS = NF90_PUT_VAR (outputID, pdf_param1_itID, pdf_param1_it)
    STATUS = NF90_PUT_VAR (outputID, pdf_param2_itID, pdf_param2_it)
    STATUS = NF90_PUT_VAR (outputID, sigma_y_itID, sigma_y_it)
    STATUS = NF90_PUT_VAR (outputID, sigmays_itID, sigma_ys_it)
    STATUS = NF90_PUT_VAR (outputID, tau_itID, tau_it)
    STATUS = NF90_PUT_VAR (outputID, nu_itID, nu_it)
    STATUS = NF90_PUT_VAR (outputID, rho_itID, rho_it)
    STATUS = NF90_PUT_VAR (outputID, y_itID, y_it)
    STATUS = NF90_PUT_VAR (outputID, acceptID, accept_vector)
    STATUS = NF90_PUT_VAR (outputID, acceptsigmayID, accept_sigma_y_vector)
    STATUS = NF90_PUT_VAR (outputID, acceptsigmaysID, accept_sigma_ys_vector)
    STATUS = NF90_PUT_VAR (outputID, accepttauID, accept_tau)
    STATUS = NF90_PUT_VAR (outputID, acceptnuID, accept_nu)
    STATUS = NF90_PUT_VAR (outputID, acceptyID, accept_y_vector)
    STATUS = NF90_PUT_VAR (outputID, HoID, H)
    STATUS = NF90_PUT_VAR (outputID, zoID, z)
    STATUS = NF90_PUT_VAR (outputID, datenumberoID, datenumber)
    STATUS = NF90_PUT_VAR (outputID, timeindexnonzerooID, timeindex_nonzero)
    STATUS = NF90_PUT_VAR (outputID, xapoID, x_ap)
    STATUS = NF90_PUT_VAR (outputID, nobsoID, nobs)
    STATUS = NF90_PUT_VAR (outputID, sitespresentoID, sitespresent)
    STATUS = NF90_PUT_VAR (outputID, deltatimeoID, deltatime)
    STATUS = NF90_PUT_VAR (outputID, distanceoID, distance)
    STATUS = NF90_PUT_VAR (outputID, TindicesoID, T_indices)
    STATUS = NF90_PUT_VAR (outputID, distributionoID, distribution)
    STATUS = NF90_PUT_VAR (outputID, DoID, D)
    STATUS = NF90_CLOSE(outputID)

  end subroutine write_data

End module MCMC_data

